Plasma-Catalysis of Nonoxidative Methane Coupling: A Dynamic Investigation of Plasma and Surface Microkinetics over Ni(111)

A heterogeneous catalytic microkinetic model is developed and implemented in a zero-dimensional (0D) plasma model for the dynamic study of methane nonoxidative coupling over Ni(111) at residence times and power densities consistent with experimental reactors. The microkinetic model is thermodynamically consistent and is parameterized based on the heats of chemisorption of surface species on Ni(111). The surface network explicitly accounts for the interactions of plasma species, namely, molecules, radicals, and vibrationally excited states, with the catalyst active sites via adsorption and Eley–Rideal reactions. The Fridman–Macheret model is used to describe the enhancement of the rate of the dissociative adsorption of vibrationally excited CH4, H2, and C2H6. In combination with a previously developed detailed kinetic scheme for nonthermal methane plasma, 0D simulation results bring insights into the complex dynamic interactions between the plasma phase and the catalyst during methane nonoxidative coupling. Differential turnover frequencies achieved by plasma-catalysis are higher than those of equivalent plasma-only and catalysis-only simulations combined; however, this performance can only be sustained momentarily. Hydrogen produced from dehydrogenation of ethane via electron collisions within the plasma is found to quickly saturate the surface and even promote the conversion of surface CH3* back to methane.

: Catalytic bed characteristics calculated based on the considered catalyst and geometrical arrangement.

Characteristic
Calculation Value The vast majority of the active sites of a catalyst are located within the pores of its pellets.
Ensuring gas phase species have access to these sites is critical for heterogeneous catalysis 1 .
Plasma contains highly reactive species, such as radicals and excited states, which can also react on the catalyst active sites, leading to synergistic effects where plasma-catalysis performs better than the equivalent plasma-only and catalysis-only cases [2][3][4][5] . However, a commonly raised issue of plasma-catalysis [6][7][8][9] is ensuring the plasma is as homogeneous as possible within the catalytic bed, in order to maximize the contact between plasma species and catalyst active sites. For the case of DBD discharges, which, by nature, operate in a non-homogeneous filamentary regime populated by streamers, it needs to be determined at which conditions the streamers are able to penetrate the pores of the catalyst. Zhang and Bogaerts 10,11 , have estimated, using Particle In Cell and Monte Carlo Collision modelling, that the Debye length is an important criterion in assessing this, as streamers can only penetrate pores with a diameter larger than this length.
The Debye length, λ D , in low temperature plasma can be calculated according to: The Debye length is estimated, for all the different gas temperatures considered, based on the steady-state results of plasma-only cases, and is compared with values from literature (Tables S2 and S3). From Table S2, it is clear that the homogeneous nature of the simulation decreases greatly the population of electrons, n e , as the power input density is 2 to 3 orders of magnitude lower than what is generally considered in pulse models [12][13][14] . Consequently, in this homogeneous model, the electron density encountered is up to 4 orders of magnitude lower than typical populations in streamers (Table S3), inevitably resulting in a Debye length that is about 1 to 2 orders of magnitude lower than that of typical DBD streamers.  3 Thermodynamic consistency of surface reaction network When building microkinetic models of surface mechanisms, it is necessary to ensure that thermodynamic consistency is respected at both enthalpic and entropic levels 17,18 . Rates of forward and backward elementary reaction steps have to be calculated under thermodynamic constraints at both individual and reaction mechanism levels. The approach applied in the current work follows that presented in 19 . For every elementary reaction j considered in the model, the following relations hold between the respective forward and backward activation energies and pre-exponential factors 17 : (1) with ∆S s,j and ∆H s,j being the reaction entropy and enthalpy of a surface process, respectively, k f 0j and k b 0j being the forward and backward pre-exponential factors, E f j and E b j the forward and backward activation energy, and R the ideal gas constant (8.314 J.K.mol −1 ).
The use of thermodynamic relationships is also employed to reduce the number of adjustable parameters, via their correlation with the chosen catalyst descriptors (the heats of chemisorption of the surface species χ i ). The enthalpies and entropies of surface reactions are estimated via thermodynamic relationships with the analogous process in the gas-phase. The thermodynamic data of surface species are explicitly obtained through correlations.
The entropy of formation of a surface species, S s,i , at a given temperature is obtained from the subtraction of the gaseous translational entropy, from the gas phase entropy of formation of S5 the corresponding gas species, S g,i , according to: with S trans−3D,i being the translational entropy of the gas phase species lost upon adsorption onto the surface.
Equation (3) implies that surface species are completely immobile on the catalyst and lose all three translational degrees of freedom upon adsorption 20,21 . The translational entropy is calculated as in 22 : with h being Planck constant (6.626 × 10 −34 m 2 .kg.s −1 ).
Analogously, the enthalpy of formation of a surface species, H s,i , at a given gas temperature is obtained by subtracting the heat of chemisorption, χ i , of the equivalent gas phase species from its respective equivalent gas phase enthalpy of formation, H g,i : The standard enthalpy (H 0 ) and entropy (S 0 ) of formation of all gas phase species considered in the model are parametrised in the form of polynomial fits in function of T 0 following the NASA chemical equilibrium code format 23 : The values of the coefficients λ i for all species are provided in Table S4 with H 0 given in (kJ.mol −1 ) and S 0 given in (J.K −1 .mol −1 ).

Kinetic model parameters calculation
The rates of elementary steps are calculated through the law of mass action with the calculation depending on the nature of the considered process. For reactions involving a gas species i, its density in the gas phase, n i , is considered, whereas for surface species i, its surface density, ξ i , is used (see Table S5). Table S5: Rates calculation of the catalytic model.

Type Process
Forward rate, Units of rate constants, Langmuir-

Hinshelwood reactions
The parameters of the kinetic model are obtained using the thermodynamic constraints defined previously and additional considerations presented below. The rate constants of any forward reaction j is calculated using the Arrhenius equation: The backward rate constant is obtained via: where K eq,j is the equilibrium constant of the reaction for which it holds: Pre-exponential factors: Pre-exponential factors k f 0,j (s −1 .P a −1 ) of species i colliding with the surface via process j, are calculated using collision theory considering an empty surface 18 . The approach is applied for all adsorption processes and Eley-Rideal reactions: The pre-exponential factor calculated as above corresponds to an upper limit, which is further corrected with a sticking coefficient, s, that represents the probability of the collision being successful 18 . For all adsorption processes a sticking coefficient of 0.1 is used while for the case of Eley-Rideal steps a value of 1 considered. The pre-exponential factors for desorption and reverse Eley-Rideal processes are obtained from equation (1), with the same sticking coefficients applied, so that thermodynamic consistency is upheld.
Pre-exponential factors for forward and backward reactions involving 2 adsorbed species, are estimated using transition state theory and entropic consistency. Specifically, the reaction entropy of each elementary step is evenly distributed around the transition state theory order of magnitude of Langmuir-Hinselwood reactions (10 11 s −1 ) to satisfy equation (1).

Activation energies and reaction enthalpies:
Dissociative adsorptions are activated as they involve the cleavage of a covalent bond. For species A that adsorbs via A + 2 * ⇌ B * +C * the forward activation energy is obtained from the unity bond index-quadratic exponential potential (UBI-QEP) method 24,25 : where ∆H ads,j is the enthalpy of adsorption estimated from: Molecular adsorptions (A+ * ⇌ A * ) do not lead to bond cleavage and consequently are assumed to be non-activated. Consequently, the forward activation energy, E ads f j , is set to 0.
For the case of surface reactions (A * +B * ⇌ C * +D * ), the activation energy of the forward step is again calculated using the UBI-QEP method according to: where the surface reaction enthalpy is calculated based on the heats of chemisorption of the participating species and the reaction enthalpy of the equivalent gas phase reaction ( The activation energies of hydrogen transfer steps, namely the Eley-Rideal processes considered in this work (A * +BH ⇌ B + AH * ), are estimated using the Polanyi-Semenov equation that was proposed by Krylov in 26 by first determining the activation energy of the forward process, where E 0 = 30 kJ.mol −1 and the enthalpy of the step ∆H ER,j is: The activation energy of the reverse process, E ER b j is obtained from Eq 2. If E ER b j is found negative, the reverse reaction is assumed to be non-activated and its activation energy is set to zero, while E ER f j is set equal to ∆H ER,j to satisfy Eq 2. Most of these processes are found to be non-activated in the direction of the abstraction of an H from the surface species by the gas one (see part 7 below).

Reactivity of vibrationally excited states:
The energy contained in vibrationally excited states of methane, ethane and hydrogen is assumed to enhance the rate of the non-reversible adsorption of these species on catalyst sites, by decreasing the activation barrier of the equivalent ground state process proportionally to the excitation energy, E ν . This mechanism is discussed in several literature works 7,27 , and is considered to be contributing to positive performances obtained with plasma-catalysis of methane upgrading 3,4 , and other processes, such as ammonia production [28][29][30] .
The approach employed here follows the Fridman-Macheret model 31 that was employed for the reactivity of vibrationally excited species in the gas phase 16 . The coefficient α is calculated using the forward and backward activation energies of the ground-state adsorption at a given temperature: This coefficient α is then used as a scaling factor to reduce the forward activation barrier of the excited state dissociative adsorption, E f (vib), such that: The reduction of the forward activation barrier leads to a significant increase of the forward rate constant of adsorption (the values for ethane and hydrogen are available in the next section).

Adsorption processes
The complete set of adsorption processes are detailed here with values at 500 K presented. Table S7: Dissociative adsorptions considered in the model along with the equivalent processes involving vibrationally excited states.

Dissociative processes
Process

Surface reactions
The surface reactions following the Langmuir-Hinshelwood formalism are detailed here with data presented at 500 K.

Eley-Rideal processes
The Eley-Rideal processes considered in the model are detailed in Table S10 with data presented at 500 K. The reactions are written in the dominant direction of hydrogen abstraction from an adsorbed species. For all the reverse processes, where a molecule (CH 4 ,C 2 H 6 , C 2 H 4 , C 2 H 2 , looses an hydrogen onto an adsorbed species, the equivalent processes involving vibrationally excited modes of the molecules were accounted for, under the assumption that the same rate as the ground-state applies.  Number densities of CH 4 and its excited states and of C 2 products, as presented in Figure 2 of the main manuscript, are provided in Figure S1 in linear scale for the different cases simulated to demonstrate more clearly the approach to steady state.

Density profiles of main radicals
The densities of the radicals with the highest populations for all cases and temperatures studied are presented in Figure S2. For all radicals, the densities are highest for plasma-only cases, followed by plasma-catalysis, while for catalysis-only values are significantly lower. The profiles of H and CH 3 radicals display much less variation with time in comparison to those of methane and C 2 species (Figure 2 in main manuscript) on account of their high reactivity and their rapid formation and consumption. The density of C 2 H 5 follows the formation of ethane, with the profile resembling that of C 2 presented in the main manuscript. At 600 and 500 K the density of radicals in the catalysis-only case is seen to increase due to the thermal activation of methane. Figure S2: Densities of the most populated radicals for all cases and temperatures studied, a) CH 3 , b) H, c) C 2 H 5

Product selectivities and methane conversion
The selectivities of all products obtained at maximum methane conversion of each temperature are given in Table S11. As the catalyst retains in cases up to 10% of total carbon (e.g. at 300 K), values reported are based on product densities according to: S CxHy = n CxHy n products For each reactor case, the selectivities towards the C 2 and C 3 species showed minimal variation for the studied temperatures. Catalysis cases were selective only towards ethane, whereas traces of ethylene were observed for plasma-catalysis cases, originating mainly from the dehydrogenation of ethane by electron collision processes. Only a small effect of temperature was observed for the plasma-only case, with conversion varying minimally, and selectivities being almost identical, between 300 to 600 K. Both catalysis cases show similar results of maximum methane conversion at 300 K, with larger deviations observed at higher temperatures. The highest conversion of the catalysis-only and plasma-catalysis cases at 400 K corresponds to the situation where methane adsorption/desorption is unbalanced and the catalyst mainly acts as a sink of methane, without emitting significant amounts of products. At 500 K the conversion in the catalysis-only scenario is higher than at 600 K as at the higher temperature adsorption/desorption cycles are faster, resulting at lower values that the densities of CH 3 * and H * stabilise at. The plasma-catalysis conversion at 500 and 600 K are similar and about 2 times higher than the equivalent plasma-only cases. These are the optimal results that are attained in these simulations, however, as discussed in the main manuscript, they are only momentary values that diminish following the adsorption of H 2 onto the catalyst. Table S11: Selectivities of C 2 and C 3 at maximum methane conversion, for all scenarios and temperatures considered. The surface species densities in catalysis-only and plasma-catalysis scenarios at 300, 400 and 600 K are presented in Figures S3, S4 and S5, respectively. Qualitatively identical trends to those at 500 K presented in the main manuscript are visible for both cases. The initial stabilisation of CH 3 * and H * occurs faster at higher temperatures, while the respective densities are higher at lower temperatures, both due to the faster adsorption/desorption of methane at higher temperature. A slight increase of H * at 600 K for the catalysis-only case is a consequence of methane thermal activation: H + CH 4 → CH 3 + H 2 accelerating and creating molecular hydrogen that dissociatively adsorbs on the catalyst. The plasma-catalysis cases all are characterised by the same divergence of CH 3 * and H * , after the initial stabilisation and the onset of ethane production. For the plasma-catalysis case, the highest coverage is observed at 400 K due to the unbalanced adsorption/desorption cycles discussed in the main manuscript. The reaction pathway analysis results on C 2 H 4 * and C 2 H * 5 , despite the important role of the species in plasma-catalysis, showed little variation with temperature ( Figures S6 and S7). At all temperatures, the population of C 2 H 4 * is rather stable, with the species having a very low net rate of production, several orders of magnitude lower than its production and consumption. Its production is largely unaffected by temperature and is dominated ethylene adsorption (C 2 H 4 + * → C 2 H 4 * ) and surface coupling CH 3 * +CH 2 * → C 2 H 4 * +H * ). C 2 H 4 * mainly hydrogenates towards C 2 H 5 * via surface H-transfer (C 2 H 4 * +CH 3 * → C 2 H 5 * +CH 2 * ). At higher temperature, due to the adsorption of H 2 , the direct hydrogenation of C 2 H 4 * into C 2 H 5 * with H * becomes significant. At maximum methane conversion, the net rate of C 2 H 5 * production overlaps with its rate of production, the latter being driven by the same processes responsible for C 2 H 4 * consumption. Only a negligible impact of the excited states of ethane C 2 H 6 (ν1, 3) and C 2 H 6 (ν2, 4) is visible at high temperature. At all temperatures, the net consumption of C 2 H 5 * is dominated by its desorption as ethane. Figure S6: Reaction pathway analysis of C 2 H 4 * . Relative contribution (%) and total integrated rates (particles.s −1 ) over temperature (K). Figure S7: Reaction pathway analysis of C 2 H 5 * . Relative contribution (%) and total integrated rates (particles.s −1 ) over temperature (K).

Hydrogen transformations in plasma-catalysis
Hydrogen transformations for plasma-catalysis at 500 K at peak simulation time t start = 0.018 s are presented in Figure S8. The dissociation of methane via electron collisions is the main source of H radicals in the gas phase. The latter are very short-lived species and mostly adsorb on the catalyst (90.89% of its total consumption). H * is formed via the dissociative adsorption of methane and molecular hydrogen H 2 . The latter is even faster than the excitation of H 2 to H 2 ν(1), H 2 ν(2) and H 2 ν(3), as visible from their minute contribution to the consumption of H 2 . Ethane is the main source of H 2 through electron collisions, confirming the propensity of the catalyst to loose activity as ethane is formed (discussed in the main manuscript). Figure S8: Reaction pathway analysis of hydrogen at 500 K and at maximum CH 4 consumption for plasma-catalysis.
14 Reaction enthalpies considered in the energy efficiency calculations The thermodynamic data used used for the energy efficiency calculations are collected in Tables S12 and S13. The enthaply of formation of all molecules at gas temperature T 0 , H T 0 f , is obtained from the NASA chemical equilibrium polynomials 23 provided in Section 3, based on which the enthalpies of reaction ∆H T 0 r are also obtained.